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ABSTRACT 

The properties of unresolved protostars and their local environment are frequently inferred from 
spectral energy distributions (SEDs) using radiative transfer modeling. In this paper, we use synthetic 
observations of realistic star formation simulations to evaluate the accuracy of properties inferred 
from fitting model SEDs to observations. We use ORION, an adaptive mesh refinement (AMR) 
three-dimensional gravito-radiation- hydrodynamics code, to simulate low-mass star formation in a 
turbulent molecular cloud including the effects of protostellar outflows. To obtain the dust temperature 
distribution and SEDs of the forming protostars, we post-process the simulations using HYPERION, 
a state-of-the-art Monte-Carlo radiative transfer code. We find that the ORION and HYPERION 
dust temperatures typically agree within a factor of two. We compare synthetic SEDs of embedded 
protostars for a range of evolutionary times, simulation resolutions, aperture sizes, and viewing angles. 
We demonstrate that complex, asymmetric gas morphology leads to a variety of classifications for 
individual objects as a function of viewing angle. We derive best-fit source parameters for each 
SED through comparison with a pre-computed grid of radiative transfer models. While the SED 
models correctly identify the evolutionary stage of the synthetic sources as embedded protostars, we 
show that the disk and stellar parameters can be very discrepant from the simulated values, which 
is expected since the disk and central source are obscured by the protostellar envelope. Parameters 
such as the stellar accretion rate, stellar mass, and disk mass show better agreement, but can still 
deviate significantly, and the agreement may in some cases be artificially good due to the limited 
range of parameters in the set of model SEDs. Lack of correlation between the model and simulation 
properties in many individual instances cautions against over-interpreting properties inferred from 
SEDs for unresolved protostellar sources. 
Subject headings: stars: formation, stars: low-mass 



1. INTRODUCTION 

Young protostars are deeply embedded in molecular 
gas and dust, which leads to significant reprocessing of 
their radiation. Since most of these sources cannot be 
directly imaged, it is necessary to infer properties of the 
protostar, accretion disk, and gas envelope using other 
means. The proliferation of infrared and millimeter in- 
struments, such as Spitzer, 2MASS, CARMA, Plateau de 
Bure, and Herschel, has contributed a wealth of multi- 
wavelength data for numerous embedded protostars over 
the last decade. Spectral energy distributions (SEDs), 
generally spanning ~ 1-1000 /im, can be constructed from 
this multiwavelength photometry. However, even with a 
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reasonably complete SED, complications such as excess 
extinction, parameter degeneracies, and fitting subtleties 
elevate the derivation of source attributes to an art form. 

SEDs are typically divided into four classes using ei- 
ther the infrared spectral slope (Greene et al. 1994) or 
the effective source ( "bolometric" ) temperature (Myers 
& Ladd 1993). These classes generally serve as a proxy 
for the amount of radiation reprocessing that occurs due 
to the embedding dense gas. Regardless of the details of 
the classification scheme, objects characterized as Class 
protostars are heavily obscured by dusty envelopes such 
that most of the radiation falls in the sub- millimeter. 
Class I protostars are less obscured and may be sur- 
rounded by a massive circumstellar accretion disk. For a 
Class II object, some of the direct stellar radiation is vis- 
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Fig. 1. — Log column density for the 0, 15, 30 and 60 kyr zooms from left to right. Star locations are indicated by the white crosses. 



ible, and often the source is considered to be a pre-main- 
sequence star, which has accreted or expelled most of its 
initial envelope and thus produces little (sub) millimeter 
emission. The remaining gas lies in a thin circumstel- 
lar accretion disk. Outflow activity may affect the SED 
shape most significantly during the Class I and Class II 
phases. For Class III objects, the disk dissipates and the 
source approaches the main sequence. 

While SED classification may be straightforward, de- 
riving actual physical properties from the SEDs is more 
complicated. For reasons discussed below, it is very dif- 
ficult to assign more than a rough evolutionary state 
to a given SED. The most common means of inferring 
properties involves picking a simple analytic model for 
the protostellar properties (e.g., mass, radius, and ac- 
cretion rate) and embedding gas morphology (e.g., disk 
mass, envelope mass, radial profile, outflow cavity width) 
and then computing SEDs with a radiative transfer code. 
The studies published by Robitaille et al. (2006, hence- 
forth R06) and Robitaille et al. (2007), greatly stream- 
lined parameter derivation by providing an extensive li- 
brary of 200,000 aperture-dependent model SEDs defined 
by 14 physical parameters (including a wide range of stel- 
lar masses and evolutionary stages) and 10 viewing an- 
gles for each physical model. The studies also included 
an efficient way to fit these models to observations and 
to derive ranges of parameter values providing a good 
fit. These models have been widely used on problems 
ranging from estimations of individual stellar properties 
(Robitaille et al. 2007) to the Milky Way star formation 
rate (Robitaille & Whitney 2010). 

Despite the breadth of parameter space, the R06 mod- 
els have limits. The morphologies are axisymmetric and 
thus do not account for gas turbulence, filamentary struc- 
ture, or outflow asymmetry. The stellar age is taken as 
one parameter input, but the models do not form an 
evolutionary sequence. (Indeed, the detailed evolution- 
ary progression of an individual protostar, including its 
formation time (McKee & Offner 2010), accretion his- 
tory (Offner & McKee 2011) and their subsequent ef- 
fects on stellar evolution (Hosokawa et al. 2011) remains 
debated!). Because of the inherent observational uncer- 
tainty in age estimations and parameter degeneracies, 
SED classes cannot be confidently mapped to stages of 
evolution (e.g., Crapsi et al. 2008). The R06 models also 
neglect gas with temperatures below 30 K, which may 
affect the fits at longer wavelengths. Finally, the models 
assume that the source is a single protostar. Although 



most low-mass stars in the field are single (Lada 2006), 
evidence suggest that protostars may have higher mul- 
tiplicity fractions (Duchene et al. 2007). In addition, 
estimating protostellar multiplicity is challenging even 
using high-resolution interferometry of local star form- 
ing regions (Offner et al. 2011a). In regions more than a 
kpc distant, where each source is potentially a cluster of 
young protostars, the multiplicity is impossible to obtain 
accurately (e.g., Kang et al. 2009). 

There are also a number of complications to parameter 
derivations that are independent of underlying model as- 
sumptions. First, there is sizable parameter degeneracy. 
This can be taken into account by considering the group 
of models with reasonable fits, i.e., "good- fit" models, 
rather than any single best-fit model. Second, the source 
viewing angle has a significant effect on the source clas- 
sification. Even including an outflow cavity inclination 
variable, it is very difficult to distinguish an older edge-on 
Class II source from a younger Class I source or a Class 
I source viewed along the outflow axis from an older, 
less embedded Class II source. This type of confusion 
can be reduced by making mass envelope estimates from 
observations of the millimeter continuum (Enoch et al. 
2009) or by detecting molecular emission associated with 
dense gas (e.g., Heiderman et al. 2010; van Kempen et al. 
2009), however, this data is not always available. 

The goal of this work is to model the SEDs of sources 
in realistic, non- magnetized simulations and explore the 
challenges of parameter estimation from SEDs based on 
idealized models (e.g., R06). We use ORION, an adap- 
tive mesh refinement (AMR) three-dimensional gravito- 
radiation- hydrodynamics code, to simulate low-mass star 
formation in a turbulent, clustered environment. The 
simulations include both the effects of radiation feed- 
back from the forming stars and protostellar outflows. In 
order to obtain model dust temperatures and synthetic 
SEDs for the sources, we post-process the simulations us- 
ing the HYPERION radiative transfer code (Robitaille 
2011). Section 2 contains an overview of our numerical 
methods. In section 3 we present the SEDs for four evo- 
lutionary times and a variety of simulation resolutions, 
viewing angles and aperture sizes. We also derive best- 
fit source parameters for each SED through comparison 
with the R06 grid of models and then compare these with 
the actual simulation values. We conclude in section 5. 

2. METHODS 
2.1. 3D Radiation- Hydrodynamic Simulations 



Synthetic Observations 
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Fig. 2. — Volume rendering of gas velocity for the 0, 15, 30 and 60 kyr zooms from left to right. 



We perform self-gravitating radiation-hydrodynamic 
calculations using the ORION Adaptive Mesh Refine- 
ment (AMR) code (Truelove et al. 1998; Klein 1999). 
Our simulations adopt the same initial conditions as 
Offner et al. (2009) and numerical procedure as Hansen 
et al. (2012), which we briefly summarize below. 

The simulations are of a turbulent box with a mean 
density of 4.46 x 10 - 20 g cm -3 , a total gas mass of 185 
M , cloud length of 0.65 pc, and an initial 3D Mach 
number of 6.6. To initialize the domain, we inject energy 
at a constant rate for two shock crossing times without 
self-gravity (e.g., following Stone et al. 1998). After this 
turbulent driving phase, gravity is turned on and the 
turbulence is allowed to naturally decay. 

Sink particles are inserted in regions of the flow that 
exceed the Jeans condition (Krumholz et al. 2004). These 
"stars" are endowed with a sub-grid stellar evolution 
model that includes accretion luminosity down to the 
stellar surface, Kelvin-Helmholz contraction, and nuclear 
burning. A second sub-grid model based upon Matzner 
& McKee (1999) governs the launching of protostellar 
winds. Our sub-grid models allows us to include all ra- 
diative heating in a way that is independent of the sim- 
ulation resolution. Since we also include the effects of 
protostellar outflows, the accretion rate from the disk 
onto the star is reduced from the case without outflows, 
and the accretion luminosity is thus smaller than com- 
puted in Offner et al. (2009) by a factor of ~ 10 (Hansen 
et al. 2012). As pointed out by Bate (2011), one de- 
fect of stellar sub-grid models is that the luminosity de- 
pends explicitly on the protostellar radius. This radius 
can vary by a factor of two depending upon the details of 
the adopted evolutionary model, including the assumed 
radiative efficiency, initial seed radius or mass, and gas 
opacities. Nonetheless, this uncertainty corresponds to a 
factor of only 2 1 / 4 (< 20%) in gas temperature (Offner 
et al. 2009). 

We summarize the key aspects of our outflow model be- 
low but refer the reader to Cunningham et al. (2011) for 
the full details of the model implementation in ORION. 
In our outflow model, the wind launching velocity is 
given by the Keplerian velocity at the stellar surface, 
vk — yj GM* I r* . For high-mass protostars these veloci- 
ties can exceed 200 kins" 1 , greatly limiting the numerical 
time step of the calculation. Cunningham et al. (2011) 
circumvent this issue by capping the outflow velocity at a 
fixed fraction of the Keplerian speed. In the calculations 
we present here, we instead economize the computational 
expense by capping the launching velocity, vk <60 km 



s _1 , and only proceeding to high resolution at discrete 
times. 

Note that in our wind model the outflow momentum 
and direction are determined by the instantaneous pro- 
tostellar mass, accretion rate, and angular momentum 
vector. Once deposited on the grid, the gas evolves hy- 
drodynamically, which leads to outflow morphology and 
asymmetry similar to that of observed outflows (Offner 
et al. 2011b). 

For the dust opacities, we adopt a dust model that as- 
sumes a standard iron abundance and treats the grains as 
composite aggregates (Semenov et al. 2003). This model 
produces an extinction ratio of R v = 3.42 for dust with 
temperature T < 120K. In most of the domain, dust 
serves as the dominant coolant. However, strong shocks 
produced by outflowing gas running into ambient mate- 
rial can result in cells with temperatures exceeding the 
dust destruction temperature of ~ 1000K. In this regime, 
we implicitly calculate the gas temperatures assuming 
atomic line cooling (Cunningham et al. 2011). 

We perform the numerical calculations as follows. 
First, we run a simulation of a forming cluster for one 
freefall time with a maximum cell resolution of 128 AU 
(two AMR levels). We then "zoom" in on four different 
time slices by restarting the calculation and allowing it to 
evolve with five additional AMR levels (4 AU resolution) . 
This is similar to the procedure in Offner et al. (2011b) 
except that here all protostars receive extra resolution. 
New grid cells are added to satisfy a Jeans condition with 
J = 0.125 (Truelove et al. 1997) and a geometric refine- 
ment criterion requiring that each star is centered within 
a block of 8 3 fine cells. We refer to the high resolution 
calculations according to the time when the zoom begins 
relative to the formation of the first star in the simula- 
tion: kyr, 15 kyr, 30 kyr and 60 kyr. Extra resolution 
is added to the first time slice beginning just prior to the 
formation of the first star particle. During each zoom, 
we evolve the calculation for ~ 8 kyr, where we allow 
structure to develop in the newest cells before proceed- 
ing to the next AMR level. Table 1 gives the protostellar 
masses and luminosities at the highest resolution during 
each of the zooms. In our analysis, the oldest star is ~ 70 
kyr, so we consider only the earliest embedded stage of 
star formation in this work. 

Figure 1 shows the column density at the completion 
of each zoom. At the final time, there are 20 stars, which 
are clustered in several different regions of the domain. 
Figure 2 shows a volume rendering of the gas velocity for 
each of the zooms. As the calculation evolves, the out- 
flows of the first forming stars extend across the domain 
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TABLE 1 

Star particle properties at the highest 
resolution for each of the zooms. the stars are 
listed in order of formation. 
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a See Figure 6. Only those with sufficient luminosity 
and those that do not have a nearby, more massive com- 
panion are assigned a unique name. 

b This source is brighter than its more massive compan- 
ion (Id) because its instanteous accretion rate is higher. 

and dominate the turbulent motions. By the final time, 
the outflows are too entangled to separate individually, 
and the simulation resembles an outflow dominated clus- 
ter like NGC1333 (Swift & Welch 2008). 

2.2. Spectral Energy Distribution Modeling 

We use HYPERION, a parallel 3D Monte-Carlo dust 
continuum radiative transfer code, to compute synthetic 
spectral energy distributions for the forming stars in the 
simulations (see Robitaille 2011 for algorithmic details). 
Since the HYPERION infrastructure is designed to be 
as generic as possible, we are able to map the ORION 
AMR cells to a tree structure that can be directly read by 
HYPERION without interpolation or information loss. 
Moreover, the sources do not need to be excised from the 
domain and computed separately. For each time slice, 
HYPERION calculates the SED for a set of apertures 
centered on each specified source as viewed through the 
box. In some cases, especially for the larger apertures, 
multiple sources may contribute to the net emission. The 
HYPERION inputs are the dust density as a function 
of position, an opacity model, source luminosities and 
source temperatures. The dust temperature is calculated 
by HYPERION based upon the input sources. For con- 
sistency, we use the same aggregate dust grain opacity 
model used in the ORION simulations and adopt the co- 
efficients calculated for dust with temperature < 120K 



(Semenov et al. 2003). 

When performing the radiative transfer calculation, we 
exclude any star particles that have not yet begun burn- 
ing Deuterium. These star particles tend to have such 
small masses that their accretion luminosity is negligi- 
ble. (These stars are shown in Figure 1 for completeness.) 
We also zero the density within 8 3 fine cells immediately 
around each star. This volume contains the accretion 
and outflow launching regions. Numerical effects reduce 
the accuracy of these cells, so that they may be over- 
dense above and below the disk plane, where the wind 
would otherwise sweep out a cavity. This is discussed in 
more detail in section 4.9 and Appendix A. 

We perform the radiative transfer calculations with 10 7 
photon packets to ensure that the shorter wavelength 
fluxes are reasonably well converged. For this number, 
additional photon packets do not significantly change the 
SED. The SEDs are computed with 200 wavelengths log- 
arithmically distributed between 0.1 /im and 5000 /im. 
For each source and viewing angle, we compute the flux 
for 20 circular apertures with radii logarithmically spaced 
between 1000 AU and 20,000 AU. Each source is viewed 
from 20 perspectives, which are regularly spaced in and 

<pK 

3. SYNTHETIC OBSERVATIONS 

3.1. Dust Temperature Comparison 

HYPERION calculates the dust temperature, T, by as- 
suming local thermodynamic equilibrium (LTE), wherein 
the dust emission balances the specific energy absorption 
rate, A: 

Attk p (T)B(T) = i, (1) 

where ftp (T) is the pre-computed Planck mean mass ab- 
sorption coefficient and B(T) is the integrated Planck 
function. The right side, which depends upon the mean 
intensity, J V1 can be expressed as a function of the energy 
per photon packet, e, cell volume, V, emission time, At, 
photon path length between events, and mass absorp- 
tion coefficient (Lucy 1999): 

f°° 1 e 

A = 47tJ n v J v dv=— — J2 1k <v ( 2 ) 

On top of this solution, we impose a minimum dust tem- 
perature of 10 K. For HYPERION, stars are the only 
heating source. Dynamical effects such as shock and vis- 
cous heating are neglected by eq. 1. Viscous heating is 
much smaller than radiative heating, where any signif- 
icant contribution falls mainly within the disks (Offner 
et al. 2009). However, shocks may significantly heat the 
gas, particularly affecting outflow gas more distant from 
the star. 

In comparison, ORION makes several different as- 
sumptions in calculating the gas temperature. First, it 
assumes that the dust and gas are perfectly collisionally 
coupled and, hence, have the same temperature and ve- 
locity. ORION assumes that dust is the primary coolant 
(except in superheated regions in which the dust has sub- 
limated). To obtain the temperature, ORION solves the 

1 Geometrically, 20 regular angular spacings in spherical coordi- 
nates correspond to the vertexes of a dodecahedron. 
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Fig. 3. — Estimated dust temperature for the first forming protostar at kyr (just after the formation of the first sink particle). Left: 
Histogram of the temperature ratios for all cells in a (1400 AU) 3 region centered on the protostar. Center: Temperature slice showing the 
ORION dust temperatures. Right: Temperature slice showing the HYPERION dust temperatures. The central 16 3 finest cells have been 
excised (this is where the HYPERION density input was set to zero). 

non-equilibrium, flux- limited, radiation diffusion equa- 
tion: 



dE 
~dt 



-V-(— VE) =K R p[47rB(T)-cE]+YS n , (3) 



where E is the radiation energy density, and ftp are 
the Rosseland and Planck mean dust opacities, A is the 
flux-limiter, p is the gas density, and S n are the stel- 
lar sources. This formulation of the diffusion equation 
omits velocity dependent terms that account for the ad- 
vection of radiative enthalpy, since these are not signif- 
icant for low- mass stars. At each timestep, the equa- 
tion is solved iteratively until the temperature converges 
in each cell. The density is determined hydrodynami- 
cally and remains constant during the iterations. The 
opacities are time-lagged and also remain constant. We 
impose Marshak boundary conditions at the edge of the 
domain, which corresponds to a 10 K radiative flux. This 
ensures that the gas temperature limits to 10 K far from 
the stellar sources. Radiative diffusion is generally a good 
approximation within dense cores where the gas is fairly 
optically thick. 

Figure 3 shows a comparison of the estimated dust 
temperatures around a young protostar. In both cases, 
the heating is mainly confined to the outflow cavity. In 
HYPERION, the difference between the core dust and 
outflow dust is due to radiative beaming. In ORION, al- 
though some beaming can occur (e.g., Cunningham et al. 
2011), the elevated outflow temperature is due primarily 
to shocks. Consequently, these cells are a factor of ~ 2 
warmer than estimated by HYPERION. Similar differ- 
ences are apparent at later times. 

3.2. The Spectral Energy Distribution Zoo 

In this section we examine the evolution of the source 
SEDs as a function of time, resolution, aperture, and 
viewing angle. 

3.2.1. Resolution 

Figure 4 shows a source observed at each zooming 
stage. To allow structure to develop at each resolution, 
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Fig. 4. — A 30 kyr source with SEDs calculated at fixed inclina- 
tion and 1000 AU aperture. The curves range from low (maxlevel 
= 3, 65 AU) to high (maxlevel=7, 4AU) resolution. Consecutive 
SEDs are ~ 1.5 kyr apart. 

there is an offset of ~ 2000 kyr between the zoom im- 
ages. 2 For consistency, we impose an inner gas radius 
of 16 AU (4 cells at the finest resolution zoom) for all 
zooms. At higher resolution, the main emission peak 
in the far-infrared shifts to shorter wavelengths. Pro- 
gressively more emission escapes in the near-infrared as 
the central disk becomes more compact and the cavity 
extends closer to the star (see Section 3.3). However, 
the trends aren't entirely monotonic with resolution. At 
long wavelengths the SED is fully converged and emis- 
sion comes from larger scales. Consequently, changes in 
the sub-millimeter portion of the SED partially result 

2 New fine cells are initialized by interpolating data from the 
parent lower resolution cells. Without additional hydrodynamic 
evolution, the data on the finest level cannot be considered higher 
resolution. 
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Fig. 5. — Source at 30 kyr with SEDs at fixed inclination (same 
as in Figure 4) and at the maximum resolution viewed with five 
different aperture sizes. 

from the time evolution between each resolution snap- 
shot. The effects of time evolution and additional reso- 
lution cannot be disentangled. 

3.2.2. Aperture 

The choice of aperture can be observationally useful 
for probing emission at specific unresolved size scales. 
Since the long wavelength emission on larger scales pri- 
marily derives from the core envelope, emission differ- 
ences between large and small apertures can be used to 
estimate the envelope gas mass. For multiple apertures, 
an average density profile can also be inferred by fitting 
how quickly the emission decreases with aperture size. 
Figure 5 illustrates the monotonically shrinking emission 
at long wavelengths as the aperture size decreases. The 
near-infrared emission also decreases slightly. This is due 
to the radiation emitted closer to the star scattering to 
larger scales. 

3.2.3. Time Evolution 

Star formation models suggest a basic progression in 
the source properties over time: the envelope mass de- 
creases, the outflow cavity widens allowing more stellar 
radiation to escape, and the disk shrinks. However, ex- 
actly how this progression occurs and specifically how 
this impacts the SED in not well known. The R06 an- 
alytic prototypes for the source, disk, and envelope pa- 
rameters encompass such changes broadly, without spec- 
ifying an exact evolutionary sequence. Our simulated 
protostars, however, do follow a consistently calculated 
evolutionary progression. Asymmetric gas morphology 
and non-monotonic evolution, which results from effects 
such as variable stellar accretion, add another level of 
complication over the R06 models. In this section we 
consider the time evolution of the SEDs in our simula- 
tions and the effect on underlying properties. 

Figure 6 shows the SEDs of the forming stars at four 
different times. The first forming star appears at all 



four times, while the second and third star appear at 
three times. In a few cases, there is a secondary or ter- 
tiary star present within the observing aperture. In these 
cases, the SEDs of the two sources are nearly identical, 
so we present only the SED for the most massive. Each 
panel shows the seds for the source viewed at four dif- 
ferent inclinations (see section 4.5 for discussion). These 
SEDs correspond to the views with the highest and low- 
est bolometric luminosities and to two views near the 
median inferred bolometric luminosities. 

Figure 6 illustrates several anticipated broad trends. 
The spectral peak shifts to shorter wavelengths with 
time. This shift is partially due to more direct and scat- 
tered emission escaping at shorter wavelengths. More 
viewing angles sample this emission with time (see sec- 
tion 3.2.4). Nearly all the sources at 60 kyr show 
some emission at shorter wavelengths, even the youngest 
sources. This is because star formation within the cloud 
is somewhat clustered, as illustrated by Figure 1, and 
the interacting outflows of the older sources impact the 
environment of the younger sources (see Figure 2). 

For the sources with close companions, we find that 
the SED appearance is not directly indicative of a sec- 
ond nearby star. Without resolving the second source 
it is thus impossible to tell which of the sources have 
companions from the SED shapes alone. Disentangling 
the effects of multiplicity on the SED is an important 
problem, but it is beyond the scope of our current study. 

Figure 7 shows the wavelength-integrated flux, i.e., lu- 
minosity, derived from the SEDs for each source and the 
range between the minimum and maximum. In some 
cases the range of luminosities spans an order of magni- 
tude. For source 7, it spans three orders of magnitude. 
In almost all cases, the input luminosity falls within the 
observationally estimated range. In the three cases where 
there is large disagreement (sources 6, 11, 13), the ob- 
served protostars are near more luminous sources. Al- 
though none of these is physically within 1000 AU of 
another source, the apertures may still overlap signifi- 
cantly with that of a brighter source, increasing the net 
emission. For example, source 6 lies ~ 1300 AU away 
from the brightest source, Id. Since it is relatively close 
to Id and dim by comparison, it is unlikely that it would 
be identified as a separate source observationally. 

Increasing the aperture size has two effects on the lumi- 
nosity ranges shown in Figure 7. First, as more envelope 
gas is included, the integrated luminosity increases. For 
large apertures, i.e., > 10,000 AU, the minimum bolo- 
metric luminosity can be 2-3 times the actual protostellar 
luminosity (even accounting for the presence of a second 
protostar). The excess emission arises from the core en- 
velope, which is radiating at ~10 K (e.g., Enoch et al. 
2008). Second, if a significant envelope remains, the IR 
flux contributes less to the total luminosity for large aper- 
tures such that the difference between the minimum and 
maximum bolometric luminosities decreases. 

3.2.4. Viewing Angle 

Previous comprehensive SED modeling has highlighted 
the influence of viewing angle on the inferred luminos- 
ity, SED shape, and apparent age (e.g., Whitney et al. 
2003). We find a similar dependence here. Figure 6 il- 
lustrates SED variability with source inclination. For 
example, one view observes source lc down the outflow 
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Fig. 6. — Sources at the four different times viewed through a 1000 AU radius aperture. The output time for each column is denoted 
in the first row. Each panel contains SEDs for four different viewing angles: the SEDs corresponding to the minimum and maximum 
bolometric luminosities and two intermediate bolometric luminosities. If there is a secondary protostar within 1000 AU only the SEDs for 
the primary source are depicted; these panels are marked with an asterisk. The inset table indicates the class of each SED. 
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Fig. 7. — Bolometric luminosities for all sources shown in Figure 
6. The input luminosity (simulation) for each model is indicated 
by the red stars. The vertical range indicates the bolometric lumi- 
nosities spanned by the different views, where the 1st intermediate 
bolometric luminosity is marked by the horizontal line. The aper- 
ture radius is 1000 AU. The luminosities of sources within 1000 AU 
of the primary are also shown (e.g., lc, Id, 3c). 

cavity, making it appear much warmer and more evolved 
than its actual age of ~ 35 kyr. In most cases, however, 
there is no clear sight line to the star since the outflow 
cavity is not entirely free of dust and gas. As a result, we 
find slightly less excess IR emission than idealized models 
for this view (e.g., Whitney et al. 2003). 

Three of the sources (numbers 1, 2 and 3) are present 
during more than one zoom. The viewing angle corre- 
sponding to the SEDs with the minimum and maximum 
bolometric luminosities shown in Figure 6 is not neces- 
sarily the same across these panels. Since the sources 
are not fixed with respect to the grid, the orientation at 
which an observer would view the source, for example 
down the outflow cavity, evolves with time. However, 
the brightest and dimmest SEDs generally correspond to 
the views down the cavity and through the edge-on disk, 
respectively. 

Figure 8 shows some typical luminosity distributions 
for 20 observed viewing angles. Younger more embedded 
sources tend to have a lower luminosity dispersion. This 
trend is apparent for source 1 in Figure 6. However, 
the picture becomes more confused when the sources are 
clustered. The environment of the younger protostars is 
contaminated by outflows and heating from older nearby 
objects. 

3.3. Protostellar Imaging 

The origin of the emission contributing to the SED 
shape is apparent in images of the constituent wave 
bands. Figure 9 shows two protostars viewed in six dif- 
ferent wavelengths. The longer wavelengths pick out the 
colder envelope dust. In the mid-infrared, the cooler, 
edge-on disk becomes visible as a narrow extinction band 
perpendicular to the outflow. Towards shorter, near- 
infrared wavelengths, scattered light highlights the out- 
flow cavity. The radiative beaming caused by the out- 
flow cavity is readily apparent. As sources become older 
and the dust warms, emission increases at shorter wave- 
lengths. The images of the younger source, which is not 
edge-on, indicate a much narrower outflow cavity. 

Figure 10 illustrates the effect of increasing the simula- 
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Fig. 8. — Histogram of bolometric luminosities for the 20 different 
viewing angles for sources 4 (left), 9 (middle) and Id (right). The 
input luminosity (simulation) for each model is indicated by the 
red vertical lines. The aperture size is 1000 AU. 

tion resolution. Overall, the higher resolution has a small 
effect on the outflow morphology. However, high density 
regions, such as the disk, become denser and more com- 
pact. 

3.4. Source Classification 

In this section, we compare the source evolutionary 
stage derived from observational metrics with the actual 
source properties. There are two main, arguably compa- 
rable, ways to assign a spectral Class to a source. The 
first involves computing the spectral slope of the SED in 
the near to mid IR (Greene et al. 1994). Young sources 
have a steeply increasing slopes (a > 0.3) that flatten 
and then decline (a < —1.6) as the source evolves and 
accretes its envelope. In this work, we use the character- 
istic spectral temperature, or "bolometric temperature" , 
to determine the spectral Class. 

The bolometric temperature, Tboi, is defined as 



T bol = 1.25 x 10 _11 (i/) K 



(4) 



(Myers & Ladd 1993) where the mean frequency is given 

by 

/ vSydv 

{v) = Js^7> (5) 

where S v is the flux density. Low bolometric temper- 
atures correspond to young, dim protostars with rela- 
tively massive, cold envelopes. The bolometric temper- 
ature gradually approaches the effective stellar surface 
temperature as the surrounding envelope is accreted or 
expelled. We separate the SEDs into early and late Class 
0, early and late Class I, and Class II (e.g., Enoch et al. 
2009; see Table 2 for definitions). The exact division of 
bolometric temperatures into classes is a matter of def- 
inition; in some cases Tboi =70 K serves as the dividing 
line between Class and Class I objects (Enoch et al. 
2009; Evans et al. 2009). Here we use 100 K as the cutoff 
between late Class and early Class I objects. 
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Fig. 9. — Two protostars in the 30 kyr time slice imaged in six 
different wavebands. The first forming star in the simulation is on 
the left. 
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Fig. 10. — Two protostars in the 30 kyr time slice (the first form- 
ing star is on the left) imaged in 1mm, where the AMR maximum 
resolution increases from 65 AU (top) to 4 AU. The images are ~ 
1.5 kyr apart. 

Figure 11 shows the inferred spectral Classes for sev- 
eral sources. The oldest (lc) shows a larger range of in- 
ferred Classes, including one view corresponding to late 
Class I. This indicates that even if a source is <30 kyr 
old, which is below the estimated Class lifetime of 40- 
100 kyr (Enoch et al. 2009; Maury et al. 2011), it may 
appear as a Class I object ~5-10% of the time. Statis- 
tically, some of these apparently older sources may be 
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offset by older sources obscured by their disks, which 
appear younger. Reassuringly, the sources in Figure 11 
appear to be early Class from most inclinations. All the 
sources have a minimum bolometric temperature around 
20 K, which is reasonable given that the minimum dust 
temperature is 10 K. 

Table 2 summarizes the number of sources observed in 
each class at each time, assuming that every view corre- 
sponds to an independent observation. In reality, none of 
these sources are older than the typical estimated Class 

lifetime, so it is surprising how many apparent Class 

1 (~8%) and Class II (~1%) sources occur. If the es- 
timated Class lifetime were ~ 70 kyr or less and the 
star formation rate were constant, then we would expect 
equal numbers of early and late Class sources when ob- 
serving the cloud at the last output time. We find that 
70% are early Class sources, suggesting that the Class 
lifetime is longer (although this argument implicitly as- 
sumes that the sources spend an equal amount of time 
in the early and late stages, which may not be the case). 

In order to diminish the role of viewing angle on the 
classification, Crapsi et al. (2008) propose using a min- 
imum envelope mass of O.1M to distinguish between 
Class I (younger) and Class II (older) sources, i.e., be- 
tween young embedded protostars (Stage I in the termi- 
nology of R06) and older disk-dominated sources (Stage 
II). We can estimate the envelope mass using 

rl 2 Q 

env = -= JTjT-T- , (6) 

where d is the cloud distance, ft mm is the flux den- 
sity at 1.1 mm observed with a 30" diameter aperture, 
ftimm = 0.0114 cm 2 g _1 is the dust opacity per gram 
of gas at 1.1 mm and Bi mm (T]j) is the Planck function 
evaluated at the dust temperature (Enoch et al. 2009). 
Applying this criterion to estimate the envelope mass for 
each source, we find that all the sources have envelopes 
exceeding 0.1 M so that their evolutionary stage would 
be correctly identified as embedded protostars. The ta- 
ble illustrates the troubling problem of associating classes 
with physical stages. Some sources that are apparently 
Class I, and hence inferred to be older than the Class 
lifetime of ~ 0.1 Myr, are in fact much younger. Since 
the length of the Class and Class I lifetimes are cal- 
culated statistically, for a large number of sources the 
effect of viewing angle will be minimized. However, the 
SEDs of individual sources may be misleading, and a fair 
fraction (~ 10%) of observed sources are inferred to be 
older than they actually are. 

4. COMPARISON WITH AN ANALYTIC MODEL GRID 

In this section, we compare the known simulated source 
properties with those inferred from fitting the R06 mod- 
els to the SEDs. Rather than using the full spectral infor- 
mation, we interpolate the fluxes at commonly observed 
wavelengths. We include wavelengths 1.25, 1.65, 2.17 /mi 
to represent Two Micron All Sky Survey (2MASS) data, 
3.6, 4.5, 5.8, 8.0 /im to represent Spitzer IRAC, 24, 70, 
160 /im to represent Spitzer MIPS, and 1.1 mm repre- 
senting data from the Bolocam continuum survey at the 
Caltech Sub- millimeter Observatory (CSO). 

The R06 models assume a dust extinction ratio of 
R v = 3.6, which is slightly higher than what we adopt in 
our simulation and radiative post-processing. However, 
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Fig. 11. — Bolometric temperatures and corresponding classifi- 
cation for three sources at 30 kyr observed in a 1000 AU radius 
aperture. The vertical bar indicates the range of values for 20 dif- 
ferent views (red), while the blue star indicates the median value. 

the R06 adopt different dust models for the disk than for 
the envelope and outflow. For the densest parts of the 
disks, R06 adopt the "disk midplane" model described 
by Whitney et al. (2003). This model has been shown to 
be a good fit for the SEDs of protoplanetary disks such 
as HH 30 (Wood et al. 2002). For the remainder of the 
dust, R06 adopt the Kim et al. (1994) model, which has 
an average particle size only slightly larger than dust in 
the diffuse ISM. Neither model includes icy grain coat- 
ings, which are expected to increase the optical depth 
by a factor of - 2 (Chakrabarti & McKee 2005). The 
former disk midplane model is more similar to the Se- 
menov et al. (2003) model we use in ORION and in the 
HYPERION post-processing. The Semenov et al. (2003) 
model predicts an opacity 7.4 times larger at 1 mm and 
~ 2 times larger at 0.1 mm than the R06 dust model. 
We discuss further in section 4.6 how opacity difference 
affects our results. 

The R06 models include luminosity from the central 
star (Siess et al. 2000) and luminosity due to accretion 
from the disk onto the protostar: 



R06 adopts / acc = 1, which assumes that all the accretion 
energy is radiated away. In contrast, ORION adopts 
/acc = 0.75 to take into account that accretion may be 
radiatively inefficient and some energy may be advected 
into the star or drive a mass outflow (Ostriker & Shu 
1995). 

In the following model comparison, we consider only 
those sources with "good" fits, which we define as \ 2 ~ 
Xbest < 3^Vdata, where 7V d ata is the number of flux data 
points considered. The model with the lowest x 2 given 
by X 2 < 30iVdata is the "best-fit" model. In a few cases 
there is only a single model that satisfies these criteria, 
but in most cases there are a number of good- fit models. 
If a second source is within 1000 AU, we only analyze the 
model parameters for the primary. Figure 12 shows two 
sources with the corresponding best-fit models overlaid. 
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TABLE 2 

Number of protostars in each category if the angle views are independent observations. 



Time Early Class Late Class Early Class I Late Class I Class II 

T ho i < 50 K a 50 K< T bo i <100 K 100 K< T bol <300 K 300 K< T bol <650 K 650 K< T bol <2800 K 



kyr 20 

15 kyr 59 1 

30 kyr 53 4 2 1 

60 kyr 233 59 22 3 3 



Model: 3008651 8 
Best fit 
X 2 = 61.517 A v = 30. 
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Fig. 12. — Observed fluxes in an 1000 AU aperture for sources 
lc (top) and 2c (bottom). Black points indicate the flux in the 
2MASS, IRAC, MIPS, and BOLOCAM bands, where the triangles 
indicate upper limits imposed by the sensitivity of the instruments. 
The black line indicates the best-fit model, while the gray lines 
illustrate models with good x 2 - The best x 2 5 dust extinction, scale, 
and best-fit model number are shown at the top. 



The younger source (lc) has only upper limits for the 
2MASS and IRAC bands. In general, we find that the 
models do not produce good fits for the most deeply em- 
bedded sources (e.g., all sources observed at kyr and 15 
kyr). This is partially due to "non-detections" in the IR, 
but also due to the source faint ness. Two specific aspects 
of the R06 models make them unsuitable for fitting very 
young sources. First, the minimum model stellar mass is 
0.1 Mq. In this case, the poor-fits are reassuring since 
the most massive star at 15 kyr is 0.05 Mq. The second 
aspect is that the dust temperatures in the R06 mod- 
els do not fall below 30 K. Since the dimmest sources 
are surrounded by 10-20 K dust, they will naturally be 
colder than expected by the models. 
Figures 13, 14, and 18 show the ratio of various model 



properties to the actual simulation values. Overall, 
the R06 models utilize more than 20 parameters (14 
of which are independent) describing the protostellar 
source, disk, outflow, and envelope. Here, we restrict 
ourselves to comparing to eight fundamental parameters: 
stellar mass, stellar radius, stellar accretion rate, source 
inclination, envelope mass, disk mass, outer disk radius, 
and inner disk radius. The root mean square of the dif- 
ference between the models and actual value are given in 
Table 3. We leave a more detailed comparison to future 
work. Also note that we compare each source for a single 
view (arbitrarily oriented with respect to the disk plane) 
and even though the source inclination is considered in 
the x 2 fit, the accuracy of the inferred parameters may 
depend upon the source orientation. 

4.1. Evolutionary Stage 

The primary physical quantity that can be derived 
from an SED is the evolutionary stage of the source. 
Since the R06 models cover a wide range of evolutionary 
stages, from embedded protostars to pre-main-sequence 
stars surrounded solely by low-mass disks, we first ex- 
amine whether the models correctly identify that all the 
sources in the simulation are in the embedded protostel- 
lar evolutionary stage. 

R06 defined three main stages of evolution based on the 
value of the stellar mass, M*, and the envelope infall rate, 
M env . The latter is in fact a direct proxy for the envelope 
density, which is defined according to the rotationally 
flattened infalling envelope models of Ulrich (1976). In 
this classification scheme, Stage I sources (or envelope- 
dominated sources) have M env /M* > 10 _6 yr _1 , Stage II 
sources (disk-dominated sources with low-density or no 
circumstellar envelopes) have M env /M* < 10 _6 yr _1 but 
M*/M* > 10 -6 yr -1 , and Stage III sources (sources with 
little circumstellar material) have M env /M* < 10 _6 yr _1 
and M*/M* < lO^yr" 1 . 

We find that all synthetic sources well fit by the 
R06 models have good-fit models with M env /M* > 6 x 
10 -5 yr -1 . Thus, the R06 models clearly indicate that 
the simulated sources are embedded Stage I sources. In 
other words, none of the models providing a good fit are 
disk-only models viewed at an extreme viewing angle. 

Note that for the purposes of the model fitting there is 
no direct luminosity associated with M env ; instead mod- 
els with high values of M env indicate the presence of a 
very dense envelope. The accretion rate associated with 
the output luminosity is M*, the rate at which mass is 
dumped from the disk onto the source (see comparison 
in section 4.4). For a realistic protostellar evolutionary 
sequence, there is likely a strong correlation between the 
infall rate and the stellar accretion rate. However, in the 
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R06 model space these two rates are independent. 

Physically, the infall rate may be many times higher 
than the stellar accretion rate if there is mass removal 
due to protostellar outflows and multiplicity (such that 
the infalling gas is distributed among two or more pro- 
tostars). For example, Hansen et al. (2012) finds a core 
efficiency factor of 1/3, which means that 2/3 of the in- 
falling core gas is expelled by outflows and thus M* is 
correspondingly lower. In Offner et al. (2010), infall rates 
in radiation-hydrodynamic simulations without outflows 
fall in the range of - 10 -6 M^r -1 - 6 x 10" 5 M^r -1 . 
These are at least a factor of 2-20 times lower than the 
infall rates suggested by the R06 best-fit models, which 
predict ^ 0.2-3 x lO _4 M yr _1 . If the cores are turbu- 
lent, as they are here, and not simply undergoing freefall 
collapse then the infall rate may be expected to be a 
factor of ~ 2 lower than suggested by the Ulrich (1976) 
model (McKee & Tan 2003). Given these very high pre- 
dicted infall rates, the stellar accretion rate must be sig- 
nificantly lower in order to avoid over-estimating proto- 
stellar luminosities. We will show that this is the case in 
section 4.4. 

4.2. Stellar Mass 

Figure 13 illustrates that inferred masses tend to be 
more accurate for the older sources, where the true 
masses lie within or very near the model range (e.g., lc- 
6). The masses tend to be overestimated for the younger 
objects, which are somewhat contaminated by dust 
heated by neighboring protostars. Since the R06 mod- 
els only include models with stellar masses > 0.1 M , 
they will tend to over-estimate the masses of sources near 
this limit or below by construction. Overall, R06 best- 
fit models achieve the best accuracy for bright sources 
in non-clustered regions. However, only three of the 13 
sources have best-fit models that disagree by more than 
a factor of 2. Given that the models span two orders of 
magnitude in source mass, the agreement between simu- 
lation and model is reasonable. 

4.3. Stellar Radius 

In contrast to the stellar mass, the R06 models sys- 
tematically overestimate the stellar radius, generally by 
a factor of ~2-3. This is understandable since the R06 
models use stellar properties interpolated from the non- 
accreting isochrones of Siess et al. (2000). In contrast 
to Palla & Stahler (1991, 1999), who agree well with the 
evolutionary model we adopt in ORION, the Siess et al. 
(2000) models do not include accretion. Instead, all pro- 
tostars are initialized with their final stellar mass and a 
stellar structure determined by hydrostatic equilibrium 
for central temperatures < 10 6 K, i.e., pre-deuterium 
burning. Consequently, the Siess et al. (2000) protostars 
begin with inflated radii > 4 R and contract towards 
the zero age main sequence. 3 For Palla & Stahler (1999), 
protostars with accretion rates of 10 ~ 5 M yr _1 termi- 
nate their accretion phase after 0.1 Myr on an L — T 
line in the H-R diagram referred to as the stellar birth- 
line. This birth-line coincides fairly closely with the Siess 
et al. (2000) H-R placement of stars with ages of 0.1 Myr 

3 In the context of the R06 models, the "initial" radius is given 
at 10 3 years, the earliest source age included in the models. 



TABLE 3 

Root Mean Square difference between the actual 
orion and r06 best-fit values for each parameter and 
the median ratio of the best-fit r06 and orion values. 



Parameter a 


RMS b 


Median(R06/ORION) 




0.4 M 


1.1 




2.8 R 


2.4 


dM*/dt 


6.8 x 10- 7 M yr- 1 


2.2 


Inclination 


18.5 deg 


0.7 


Menv 


0.9 M 


3.2 


A^env,4840AU 


57.0 M 


6.2 


M d 


0.02 M 


11.2 


Rd 


174.1 AU 


0.3 


Rd,in 


27.6 AU 


0.3 



a All parameters are compared assuming a 1000 AU aperture 
except M enVj 4840AU, which is calculated using a 4840 AU aper- 
ture. 

b The RMS is defined as RMS = [Oi ~^ 06i]2 ^j . 

and masses M < 1M (e.g., Da Rio et al. 2010). Conse- 
quently, the two sets of models give similar values for the 
stellar radii as the stars approach 0.1 Myr. For stars with 
ages < 0.05 Myr, the Siess et al. (2000) models overesti- 
mate the radii relative to our ORION model by a factor 
of ~2. Palla & Stahler (1991) find that assuming lower 
accretion rates, e.g. 1O _6 M yr -1 , produces smaller pre- 
birth-line radii. Thus, for protostars with low accretion 
rates and young ages, the two models will be even more 
discrepant. Some of the difference in model radii is offset 
by the ORION assumption that the accretion energy is 
not perfectly radiated away, i.e., / acc = 0.75. This re- 
sults in an effective stellar radius of (l// a cc)-R*> which is 
somewhat closer to the estimated values of Siess et al. 
(2000) 

Despite clear model differences, determining the ab- 
solute correctness of the stellar evolutionary models is 
not straight-forward. The initial protostellar radius, 
its subsequent evolution, and how it depends upon ac- 
cretion and radiative efficiency, remains hotly debated 
(Hosokawa et al. 2011; Baraffe et al. 2009). Even without 
considering these effects, various calculations often adopt 
initial radii of 2.0R o (Stahler 1988) to 3.5R© (Palla & 
Stahler 1991). In cases where accretion is considered, ac- 
cretion rates comparable to ~ lO _6 M yr _1 , which are 
appropriate for low-mass stars, tend to have early radii 
of - 1.5i? (Palla & Stahler 1991; Tan & McKee 2004; 
Hosokawa & Omukai 2009). Given the differences be- 
tween stellar models, it is encouraging that many of the 
protostellar masses inferred are nonetheless quite close 
to the simulated values (see section 4.2). 

4.4. Stellar Accretion Rate 

The R06 models achieve fairly good agreement with 
M*. With the exception of two cases, the models either 
bracket the actual value or come very close to it. Since 
the good-fit models may span four orders of magnitude 
in the accretion rate, we plot the fourth root to reduce 
the error bar extent in the plot. 

In the R06 models, the range of model accretion 
rates decreases monotonically with source age. If proto- 
stars experience episodic accretion (Hartmann & Kenyon 
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Fig. 13. — Ratio of the inferred best model values to the actual 
simulation value for each of the sources with good fits, where M* 
is the stellar mass, R* is the stellar radius, (dM* /(it) 1 / 4 is the 
fourth root of the accretion rate onto the star, and the Inclination 
refers to the tilt with respect to the line of sight. The dotted line 
indicates where the models determine a value identical to the true 
value in the simulations. Points with no error bars have only one 
best fitting model parameter. 
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Fig. 14. — Ratio of the inferred best model values to the actual 
simulation value for each of the sources with good fits, where 
is the disk mass, Rd is the stellar radius, and M en v is the envelope 
mass. The dotted line indicates where the simulation and models 
are equivalent. Points with no error bars have only one best fitting 
model parameter. Boxed points indicate that the comparisons use 
the simulation upper limits for the disk mass or radius. 

1996; Zhu et al. 2009; Vorobyov & Basil 2010; Offner et al. 
2011b), then protostellar accretion rates may jump be- 
tween 1O -7 M yr _1 and 1O -4 M yr _1 , extremes which 
may not be spanned by the models and thus result in 
age or mass errors. However, we find that, at least for 
these early times, variability in the instantaneous accre- 
tion rate is less than an order of magnitude. Conse- 
quently, our simulated values fall within the model ranges 
of R06. The extreme accretion events known as FU Ori 
outbursts are believed to result from a combination of 
thermal and magnetic instability in the inner accretion 
disk (Zhu et al. 2009). We do not treat these effects in 
these simulations, and thus all accretion variability is due 
to gravitational instability within the disk or turbulent 
variations of the core envelope. 
The good agreement between the stellar accretion rate 



determined from the SED models and the accretion rate 
in the simulation is surprising, nonetheless, since the HY- 
PERION source modeling assumes a single input spec- 
trum with a bolometric luminosity equal to the sum of 
the stellar and accretion luminosities. Note that this sum 
is ultimately what sets the range of possible good-fit val- 
ues for the protostellar masses and accretion rates shown 
in Figure 13. However, since there are many different 
values for the stellar and accretion luminosities that give 
the same bolometric luminosity, it should not be possible 
to disentangle the relative contributions from these from 
SED modeling. Instead, it is likely that the accretion rate 
is only constrained compared to the original range in the 
model grid due to other parameter constraints. In par- 
ticular, due to the parameter space sampling in the R06 
grid, models with higher envelope infall rates have higher 
stellar accretion rates, and conversely, models with lower 
envelope infall rates have lower stellar accretion rates 
(see Figure 1 of Robitaille 2008). Thus, constraining the 
envelope infall rate (or density), automatically restricts 
the stellar accretion rate to a range that is smaller than 
that of the whole set of models. Similarly, constraining 
the total luminosity also automatically rules out models 
with stellar accretion rates that are too high. 

Figure 15 shows the distribution of stellar accretion 
rates for all models in the R06 set, for the models pro- 
viding a good fit, and for all the models in R06 that have 
envelope infall rates and bolometric luminosities within 
±10% of the values of the good fits. 4 The fi sure demon- 
strates that while the range of accretion rates for the 
good fits is much narrower than the range in the entire 
set of models, this range is primarily set by constraints 
on other parameters. It is therefore likely that M* only 
appears to be well fit due to constraints on other param- 
eters that can be more directly determined, such as the 
envelope infall rate and the bolometric luminosity. These 
two parameters supply the main constraints on the stel- 
lar accretion rate either by setting an upper bound (the 
luminosity) or by restricting the range for given enve- 
lope properties (the infall rate). However, we note that 
the agreement does appear to indicate that the accretion 
rates assumed in the R06 models span a realistic range. 

4.5. Inclination 

Inferring the correct source inclination is important 
for accurately determining other source properties. If 
a source is observed edge-on, but not identified by the 
models at this orientation, then it will appear erroneously 
younger. We use the simulated angular momentum of the 
protostar to determine the inclination with respect to the 
line of sight, where the net angular momentum vector is 
assumed to point perpendicular to the disk plane. Figure 
13 illustrates that the most discrepancy occurs where the 
inclination is misfit, e.g., stars 8 and 9. These sources are 
inferred to have higher stellar masses and lower accretion 
rates than they actually do. At intermediate inclinations 
between pole-on and edge-on the inclination parameter is 
fairly degenerate. As a result, the inclination constraints 
are ultimately not very tight. 

4 Since the infall rates and bolometric luminosites have been 
varied independently of the other fit parameters, these models are 
not necessarily close, i.e., within 10% of being good-fit models, to 
the data. 
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4.6. Envelope Mass 

For each R06 model, the envelope radius, i? e nv, is ran- 
domly drawn from a uniform logarithmic distribution of 
values between Rq x 4 and i?o/4, where Rq is approxi- 
mately the radius at which the optically thin radiative 
equilibrium temperature declines to 30 K: 
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These values are further constrained to lie between 10 3 
AU and 10 5 AU. Inside this radius, the density pro- 
file goes as p(r) ex r -3 / 2 with some slight modifica- 
tion due to rotational flattening. The envelope mass 
is then computed by integrating the mass inside R en v 
The ambient gas outside the envelope is assumed to be 
a constant density randomly drawn from the range 50 
cm" 3 < n H2 /(M*/M ) < 100 cm" 3 . Thus, the envelope 
mass is coupled to both the density distribution and stel- 
lar properties. 

In the simulations, however, the gas morphology is 
independent of the stellar properties and is not nicely 
spherically distributed. The mass asymmetry means that 
there is no clear cut way to define a physically meaning- 
ful envelope mass. Instead, we follow the R06 model 
convention and calculate the total mass inside a spher- 
ical volume with the radius R env suggested by the best 
fit model. This circumvents the need to define a density 
minimum or an effective radius for the simulations. 

Figure 14 shows that the model envelope masses run 
high relative to the simulated envelope masses. In very 
few cases do the best fit models encompass the true value. 
This is primarily a result of the different opacity models 
adopted by ORION and R06. The envelope mass is de- 
termined from the long wavelength emission. At 1mm, 
the Semenov et al. (2003) opacity is 7.4 times higher than 
R06 dust opacity; it is ~ 2 times higher at 0.1 mm. Since 
the envelope mass is inversely proportional to the opacity 
(see eq. 4.6), a given flux will be interpreted by the R06 
models to correspond to a larger mass envelope. Figure 
16 shows that the discrepancy is less for a 1000 AU aper- 
ture than for the larger apertures, which include more of 
the cold, extended envelope. 

However, in all cases there is a fairly large scatter in 
the best-fit models. Since the simulated cores are embed- 
ded within a cloud and somewhat clustered, it is likely 
that some of the scatter is due to difficulty separating en- 
velope emission from foreground emission, and thus cor- 
rectly determining the line-of-sight extinction. The value 
of the inferred envelope radius is also important since 
most of the envelope mass is found a large radii. Since 
the R06 models only include envelope dust with tem- 
peratures > 30 K, they preferentially neglect the coldest 
envelope mass. This can contribute to error in the 1 mm 
flux data point, which influences the determination of the 
envelope mass. Despite this innaccurary, the models do 
correctly identify the SEDs as belonging to embedded 
protostars and do not mistake them for more evolved 
disk-dominated sources. 

For a given observed SED there will always be some 
underlying uncertainty in the grain composition and size 
distribution, which depends upon the local temperature, 
metallicity, and source age. However, our comparison 
suggests that without icy coatings and inclusion of cold 
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Fig. 15. — The top and bottom panels show, for two different 
sources, the distribution of stellar accretion rates, M*, in the entire 
R06 set of models (black), the models that provide a good fit to 
the SED (green), and the models in R06 that have envelope infall 
rates and bolometric luminosities within ±10% of the values for 
each good fit. The vertical dashed lines indicate the best-fit model. 



envelope gas, the R06 models may systematically over- 
estimate observed envelope masses. 

4.7. Disk Mass 

We follow the same procedure as Offner et al. (2010) to 
define what constitutes a disk in the simulations. First, 
we identify all connected cells around the protostar with 
densities > 10 -15 g cm -3 . This threshold effectively se- 
lects for gas comprising the disk, since the density of the 
gas in the envelope falls off rapidly. We then estimate 
the total angular momentum vector of the gas and rotate 
the coordinate frame so that the net angular momentum 
vector is parallel to the z-axis. Finally, we restrict the 
vertical disk height in the z direction to ±5 cells from 
the disk mid-plane. This is to eliminate confusion with 
streams of gas feeding the disk. The disk mass is then 
the integrated gas mass in these cells. 

Due to the inclusion of outflows, we find that the disks 
are less orderly than in the simulations analyzed in Offner 
et al. (2010). In a few cases, i.e., where the protostellar 
mass is < 0.1 M , there is no gas exceeding the density 
threshold. For these, we assume that the disk is too 
small to be resolved by the simulation and set the disk 
radius, to 4 fine cells (16 AU), the radius of the 
outflow launching region, and the disk mass to the mass 
enclosed in a spherical volume with this radius. These 
values thus represent upper limits on the disk mass and 
radius. 
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Fig. 16. — Ratio of the inferred best-fit envelope mass to the 
actual envelope mass for each of the sources with good fits as a 
function of aperture size. The dotted line indicates where the sim- 
ulation and models are equivalent. The dashed line indicates the 
ratio of the ORION opacity to the R06 model opacity (7.4) at 1 
mm. All inferred envelope masses are shown (not only those cen- 
tered on the most massive protostar in a multiple system). Each 
source has a unique symbol. 

Figure 14 shows that the model disk masses are not 
well correlated with the simulated disk masses. The mod- 
els tend to overestimate the disk mass, and in some cases 
the discrepancy is a factor of 10. Unlike other parame- 
ters, such as the stellar mass, the actual disk mass does 
not fall within the range of good-fit models in nearly all 
cases. This is mainly because the simulated disks are 
more complex and less dense at the mid-plane, where 
the grid resolution does not resolve the scale height. The 
mass estimates are most overestimated in cases where 
the R06 model underestimate the size of the inner hole 
(see section 4.9), and thus overpredict the mass in the 
inner disk region. Since the simulated disks are more ex- 
tended and less dense, it is also challenging for the models 
to accurately distinguish between emission from the disk 
and emission from the envelope, which contributes to the 
scatter in the agreement. 

4.8. Disk Radius 

In the R06 models, the disks are assumed to be per- 
fectly symmetric with a sharp edge, so there is no diffi- 
culty in defining the disk radius. In the simulations, the 
disks are often asymmetric and occasionally have spiral 
structure that winds out to larger distances. For the 
simulations, we define the radius as the average radial 
distance of all the cells in the disk. This minimizes the 
effect of asymmetry on the radius estimate. 

Figure 14 shows that the model disk radii are generally 
much smaller than the simulated disk radii. The best- 
fit model disks have a median radius of ~32 AU with a 
minimum value of 5.7 AU, while the simulations have a 
median of ~136 AU and minimum (upper limit) of 16 
AU. This is to be expected since the simulations do not 



resolve the disk scale height and thus likely over estimate 
the disk extent. Since the disk mass is distributed over a 
larger volume and the gas is less dense than expected by 
the R06 models, some of the gas likely resembles envelope 
gas for the purpose of model fitting. However, the large 
range of radii inferred by the good-fit models suggests 
that this parameter is not well constrained by the SED 
in general. The sources where the model agrees within 
a factor of two are all comparisons assuming the upper 
limit of 16 AU for the simulated disk radius. In princi- 
ple, the disk properties should be sensitive to the viewing 
angle, however, we don't find that sources with better fit 
inclinations have necessarily more accurate disk proper- 
ties. The accuracy of the inferred disk properties also 
likely depends on evolutionary stage. For more evolved 
sources, where the envelope mass is small, the R06 may 
determine more accurate disk masses and radii. 

4.9. Disk Inner Radius 

We choose to ignore the gas inside the accretion region, 
which corresponds to the 8 3 cells centered on each star. 
This effectively enforces that all accretion disks have an 
inner radius, or hole, of ~ 16 AU. Since gas is removed 
from these inner cells during the accretion process at each 
timestep the dynamics in these cells is predominantly 
numerical. The outflow is launched just outside these 
cells so that gas remains close to the star when otherwise 
it might be evacuated by the outflow. To account for this, 
we set the density to zero in these cells when observing 
the protostars with HYPERION. 

For Class sources, which are heavily embedded, the 
details of the inner region have minimal effect on the 
SED characteristics. Figure 17 shows three SEDs cal- 
culated with and without the accretion region gas. In- 
cluding the accretion gas produces an increased flux from 
~ 8 — 20 /im, i.e., increased photon scattering off of the 
inner warm dust. In the early Class case (e.g., the 
green curve in Figure 17), the dust is optically thick to 
< 20/im emission and the inner region has negligible in- 
fluence. The SEDs of late Class and I sources are more 
sensitive to this emission. In Figure 17 the bolometric 
temperatures change by 3.6% (late Class 0) and 7.8% 
(late Class I), while the early Class bolometric tem- 
perature only increases by 0.2%. However, the difference 
turns out to have little effect on the properties inferred 
in the model comparison, which uses only discrete wave- 
length points from the SED. The largest flux difference 
falls between the observed 8/im and 24/im data points. 
Consequently, the inferred model properties are largely 
insensitive to the gas in the accretion region. (See Ap- 
pendix A for additional discussion.) 

We can verify that this is the case by comparing the 
inner disk radius inferred by the best fit models with the 
actual inner disk radius of 16 AU. The R06 models in- 
clude disks with inner holes as large as 100 AU in order 
to account for the possibility of a close, unresolved stel- 
lar companion that is evacuating the inner disk. Figure 
18 shows the ratio of the inferred and actual inner disk 
radius. The large range of good- fit radii and large dis- 
crepancies of the best fit model suggest that the SED 
shape is not particularly sensitive to the inner hole size. 
At least for Class objects, the accuracy is also insensi- 
tive to the magnitude of the bolometric temperature. 
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Fig. 17. — Three SEDs calculated with (dashed) and without 
(solid) the gas inside the accretion region. The three classifications 
are early Class (green), late Class (red) and late Class I (blue). 



T b0 , (K) 



10 20 30 40 50 60 



1c 2b 2c 3c 4 5 6 



9 10 11 



Fig. 18. — Ratio of the inferred best model inner disk radius to 
the actual simulation value for each of the sources with good fits. 
The dotted line indicates where the models determine a value iden- 
tical to the true value in the simulations. Points with no error bars 
have only one best fitting model parameter. The color indicates 
the bolometric temperature. 

5. CONCLUSIONS 

We present results from turbulent, gravito-radiation- 
hydro dynamics AMR simulations of forming stars. By 
including the effects of protostellar outflows and radia- 
tive heating, we are able to follow the evolution of the 
protostars self-consistently using realistic physics. We 
neglect magnetic fields, which we leave to be addressed 
in future work. 

We use the Monte-Carlo radiative transfer code HY- 
PERION to compute the synthetic spectral energy distri- 
butions of the sources for a variety of apertures, viewing 
angles, and resolutions. We find that the HYPERION 
and ORION calculated dust temperatures are most dis- 



crepant in shock heated gas such as within the outflow 
cavity, but generally agree within a factor of two for most 
cells. 

Synthetic images of the forming protostars in the in- 
frared and (sub) millimeter show clearly recognizable ob- 
servational features. Emission from ~l-30 /im traces out 
the outflow cavity. In the sub- millimeter, an edge-on 
dense disk appears as a dark band obscuring the proto- 
star. 

Like previous work, we find that the source SED and 
its inferred properties are very sensitive to the viewing 
angle. The inferred bolometric luminosity variation can 
exceed a factor of 5 depending on the view. This corre- 
sponds to a typical standard deviation in log luminosity 
of 0.2. However, in most cases, the inferred range in- 
cludes the true luminosity. The variety of SEDs also 
corresponds to a large range of bolometric temperatures, 
which can often span two spectral classifications. This 
underscores the difficulty of associating classes with spe- 
cific evolutionary states based on SED characteristics 
alone. Since star formation is clustered, we find that 
the later forming stars suffer environmental contamina- 
tion from the neighbors, increasing the uncertainties in 
their classification and inferred properties. 

We next assess the accuracy of source properties in- 
ferred from analytic models. We use the model grid of 
Robitaille et al. (2006) to derive best-fit parameters for 
each of the sources and discuss the accuracy for a num- 
ber of fundamental parameters. No good fits occur for 
the dimmest, most embedded objects, i.e., those with 
ages < 20 kyr. The sources well fit by the R06 mod- 
els are correctly identified as embedded protostars and 
not mistakenly fit with more-evolved disk-only models. 
Using the R06 models, we thus come to the same con- 
clusion as when assessing the observed envelope mass: 
all the protostars are embedded Stage I objects. Since 
the R06 models do not include stars with final masses 
below 0.1 Mq, which is the critical envelope mass, there 
is small advantage to using the R06 models purely to 
determine the evolutionary stage in lieu of simply esti- 
mating the envelope mass. However, the R06 models 
offer many additional, potentially enlightening, parame- 
ter constraints. Further work extending the R06 model 
grid to lower masses may facilitate the evolutionary char- 
acterization of lower mass, young sources. 

The model parameters determined from the good SED 
fits for the older sources exhibit varying amounts of ac- 
curacy compared to the true simulation parameters. We 
find that the stellar accretion rate spans the true value 
in most cases. The models are consistent with the true 
stellar mass for the older and/or more isolated sources. 
This is partially because the R06 models have a mini- 
mum mass of 0.1 Mq, and thus, will not well fit young 
low-mass sources by construction. In several cases, where 
the source environment is affected (via heating or outflow 
activity) by a more luminous neighbor, the stellar mass 
estimations are quite discrepant. The R06 models also 
use the Siess et al. (2000) stellar evolution models and 
thus systematically overestimate the source radii at early 
times. This may also introduce error into the parameter 
estimates. The inferred source inclination exhibits de- 
cent agreement with the actual inclination, although this 
parameter is not strongly constrained by the fits. The 
models systematically overestimate the envelope mass, 
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a trend that might be corrected in future models that 
include icy grains and do not exclude the coldest enve- 
lope gas. We find that the disk mass and radius are not 
well fit by the models. While we expect that including 
magnetic effects and even better resolution would lead to 
smaller disks, ideally the breadth of the models should 
encompass the full disk parameter space. Some numer- 
ical studies have found that magnetic effects may com- 
pletely suppress the formation of disks (Galli et al. 2006; 
Price & Bate 2007; Hennebelle & Ciardi 2009), which 
would also be challenging for the R06 models since they 
assume that all young sources have a disk. 

Overall, our study highlights a number of problems 
with inferring source, disk, and envelope properties from 
observed SEDs of protostellar sources. Consequently, the 
classification and inferred properties of individual sources 
should be accepted with caution. Properties of isolated 
protostars are likely to be more accurate than protostars 
in clustered regions. As summarized in Table 3, proper- 
ties such as the inferred stellar mass may be inaccurate 
by more than a factor of two and the envelope mass may 
be overestimated on average by a factor of four. The 
comparison is performed assuming that the sources are 
unresolved, so additional information from direct imag- 
ing might improve the accuracy of inferred parameters. 

Some of our comparisons underscore uncertainty in 
fundamental star forming properties, such as dust prop- 
erties and the dependence of stellar evolution on accre- 
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Fig. Al. — Distribution of optical depths at 8/mi for each of the viewing angles for the sources at 60 kyr. The vertical lines indicate the 
approximate divisions between classes (see Table 2). 
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APPENDIX 

SENSITIVITY OF THE RESULTS TO THE INNER DISK RESOLUTION 

With or without the inclusion of the accretion region gas, the inner disk properties are not well resolved even with 
4 AU cell resolution. Observationally, the warm inner disk gas contributes SED emission at < 8 fim wavelengths 
(Robitaille et al. 2006). In observed SEDs, this emission can be used to directly infer disk properties (e.g., disk mass, 
flaring, gap size) of non-embedded protostars (e.g., Espaillat et al. 2011). In instances where the protostar is deeply 
embedded, namely during the Class phase, radiation from gas close to the star is significantly reprocessed by the 
cold, envelope gas. As we have demonstrated in section 4.9 the properties inferred by the R06 models are relatively 
insensitive to the details of the inner disk in most cases. 

We can estimate the importance of the inner regions more quantitatively by estimating the 8/im optical depth for 
different protostellar views. In cases where the optical depth t% < 1, our models will not accurately represent the SED 
at short wavelengths. For ~ a few, the gas is optically thick to the emission from the inner regions, and the SEDs 
will not be significantly effected by the details of the numerical simulations on these size scales. 

Figure Al shows the distribution of optical depths at 8 /am for all the different views of the protostars at the final 
output time. Only 27 views have optical depths < 1. This approximately corresponds to the number of Class I and II 
sources that we find. Since these objects are viewed on sight-lines down the outflow cavity, we would a priori expect 
a low optical depth. The large majority of sources are Class 0. Based upon these optical depths, we expect that our 
simulated early Class and most late Class SEDs accurately represent real SEDs. The simulated SEDs of later 
classes will only accurately represent protostellar disks with large inner holes. Our results in section 4, which are 
obtained for Class sources, should thus be generally applicable and relevant to observations. 



